Moments of spectral functions: Monte Carlo evaluation and verification 
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The subject of the present study is the Monte Carlo path-integral evaluation of the moments of 
spectral functions. Such moments can be computed by formal differentiation of certain estimating 
functionals that are infinitely-differentiable against time whenever the potential function is arbi- 
trarily smooth. Here, I demonstrate that the numerical differentiation of the estimating functionals 
can be more successfully implemented by means of pseudospectral methods (e.g., exact differen- 
tiation of a Chebyshev polynomial interpolant), which utilize information from the entire interval 
(— /3H/2, ph/2). The algorithmic detail that leads to robust numerical approximations is the fact 
that the path integral action and not the actual estimating functional are interpolated. Although the 
resulting approximation to the estimating functional is non-linear, the derivatives can be computed 
from it in a fast and stable way by contour integration in the complex plane, with the help of the 
Cauchy integral formula (e.g., by Lyness' method). An interesting aspect of the present development 
is that Hamburger's conditions for a finite sequence of numbers to be a moment sequence provide 
the necessary and sufficient criteria for the computed data to be compatible with the existence of 
an inversion algorithm. Finally, the issue of appearance of the sign problem in the computation of 
moments, albeit in a milder form than for other quantities, is addressed. 
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I. INTRODUCTION 

Perhaps one of the most outstanding insufficiencies 
of the path integral formulation of quantum mechan- 
ics is that it does not lead directly to efficient algo- 
rithms for the computation of dynamical information. 
In contrast, statistical quantities or imaginary-time data 
are relatively easy to evaluate and the mechanisms to 
do so are well understoodi^*^ In principle, it is possi- 
ble to relate the quantum correlation functions in real 
time to their imaginary-time counterparts by analyti- 
cal continuation. 4,5 However, such attempts lead to in- 
verse problems which, albeit uniquely determined, are 
fundamentally ill-posed. An example is represented by 
the real inverse Laplace transform^ a technique that 
has been extensively utilized as a link between the real 
and imaginary-time worlds. It requires the resolution 
of a linear integral equation that becomes extremely ill- 
conditioned upon discretization. The treatment of such 
problems is the domain of regularization theory, which 
attempts to stabilize the resulting equations by control- 
ling certain properties of their solution. 

Notwithstanding earlier attempts that were more or 
less in tone with Tikhonov's least-square approachj§iL£ 
the pioneering research of Gubernatis, Jarrell, Silver, and 
Sivia^ has spearheaded the application of methods of 
Bayesian statistical inference with entropic priors 5 ^ as a 
regularization technique for the inversion of the Laplace 
transform. To give a few examples, although somewhat 
restricted to the direct research interests of the present 
author, such methods have been successfully applied for 
the computation of various physical properties such as 
spectra*i2iiAii2ii£ quantum rates of reaction? 1 ^ 5 - and dif- 
fusion constantsiiS. When we say "successfully," we take 



into account the fact that, in most cases, there are vir- 
tually no computationally feasible alternatives: the tech- 
niques based on imaginary-time data are amenable to 
direct Monte Carlo path integral treatment and exhibit 
little degradation of their stability with the increase in 
the physical dimensionality. This is so because the sta- 
bility is related to the properties of the spectral function 
(a one-dimensional probability distribution) and not to 
the dimensionality of the physical system. 

Nevertheless, as Jarrell and Gubernatis point out? 5 " "to 
solve an ill-posed problem, nothing beats good data." 
The present paper does not attempt to improve on pre- 
vious results regarding the stabilization of the inverse 
problems. Rather, its purpose is to provide a means 
to obtain high-quality input data for the reconstruction 
of various autocorrelation functions for physical systems 
in the continuum space. A recent study of the present 
author 1 ! suggests that one way to achieve better results 
is to break the inverse Laplace problem into two sepa- 
rate steps: computation of moments by differencing an 
estimating functional (which is related to the imaginary- 
time correlation function) followed by resolution of the 
ensuing symmetric Hamburger moment problem. Both 
steps are exponentially unstable, albeit to a lesser degree 
than the original problem. Very likely, their combined 
effect is a problem that is as ill-conditioned as the origi- 
nal one. Quite clearly, we cannot create new information 
in a stable and consistent manner just by manipulating 
the data. Is there any gain in the new approach? 

The reason we shy away from exponentially unstable 
algorithms is that they require the utilization of expo- 
nentially fast algorithms (polynomial in the number of 
digits) for the computation of the input data. Unfortu- 
nately, most of the algorithms we posses converge poly- 
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normally in the best circumstances and there are only a 
few instances of exponentially fast algorithms. Some of 
the more interesting examples are related to the execu- 
tion time of elementary functions on a classical computer, 
in arbitrary precision. In fact, Brentji^ has shown that 
the evaluation of elementary functions can be performed 
in time proportional to 0(n log(ro) 2 log log(n)) with re- 
spect to the number of digits n. Therefore, as Ref. ar- 
gues, differencing an arbitrarily smooth estimating func- 
tional is a computationally feasible approach whenever 
empirical potentials are utilized. The present author has 
obtained excellent results in unpublished tests that have 
employed Amber force fields 1 ^ and Bailey's arbitrary pre- 
cision package MPFUN903 

However, such a technique cannot be applied for po- 
tentials that are the result of a computation performed 
in polynomial time. In addition, there is some unease 
related to the mere requirement of arbitrary precision. 
In most circumstances, we are interested in learning the 
properties of the spectral function in the low frequency 
region (for some transport phenomena, we are interested 
in the value at the origin of the spectral function asso- 
ciated with the flux- flux 14, 15 or velocity- velocity^ corre- 
lation functions). Due to quantum and thermal smooth- 
ing, the low-frequency portion of the spectral function 
is largely insensitive to the precision with which the po- 
tential is known. Clearly, whether we use an otherwise 
smooth potential with 7 (single precision) or 15 (double 
precision) significant digits, we do not expect the value 
of a diffusion coefficient to change dramatically and this 
expectation is justified in many cases by results of per- 
turbation theory. It follows that the instability of the 
differencing step is a property of the algorithm and not 
necessarily an inner characteristic of the problem. 

In Section IV, we show that the instability associated 
with the differencing step can be removed by interpola- 
tion of the action. That is, the path-integral action, re- 
garded as a function of the imaginary time on the interval 
(— (3H/2, (3h/2), is replaced by a smooth interpolant con- 
structed by means of trigonometric or Chebyshev poly- 
nomials. For infinitely differentiable potentials, the in- 
terpolant converges faster than any polynomial and it 
rapidly feels the discontinuities due to the finite preci- 
sion in the computation of the action. However, if it 
is known that the potential is smooth, one can utilize 
a low-degree interpolant only. By the arguments in the 
preceding section, the properties of the spectral functions 
in the low-frequency region and, therefore, the values of 
the low-order moments are not sensitive to the errors in 
the action. Despite the fact that the interpolated action 
may not necessarily come from a perturbed potential, the 
numerical results of Section V suggest that the computed 
Monte Carlo data still represent a sequence of moments, 
even for low-order interpolants. 

We commence this paper, however, with Section II, 
where we present a short review of the moment prob- 
lem and discuss its relevance for the reconstruction of 
spectral functions. In Section III, benefitting from the 



existence of certain mathematical results concerning the 
Hamburger moment problem, we give necessary and suf- 
ficient criteria for inversion algorithms to exist. These 
criteria can be utilized as an a posteriori verification tool 
for the computed data and they reveal the exponential 
extent of the instability of the moment problem. Nev- 
ertheless, this instability has to be weighted against the 
fact that the information furnished even by a few tens of 
moments is, in general, quite substantial. In Section VI, 
we review the main findings of the present work and enu- 
merate several research issues left outstanding. 



II. THE INVERSE PROBLEM AND THE 
POSITIVITY REQUIREMENT FOR THE 
SPECTRAL FUNCTION 

The dynamical quantity we want to evaluate is a quan- 
tum correlation function of the typa^i 

Co At) = tr ^ e -W^+it/n)H t e - W 2-x-it/n)H ^ 

(1) 

The operator H stands for the Hamiltonian of the sys- 
tem, a self-adjoint and bounded from below operator, 
whereas t £ R and (3 — l/(fc,gT) > are the real time 
and the inverse temperature, respectively. denotes 
the adjoint of the operator O. The parameter A may 
take values in the interval (—(3/2, (3/2). The values of the 
correlation functions for A = ±(3/2 can be recovered as 
the corresponding limits, by continuity arguments. Many 
quantities of physical interest are related to the quantum 
correlation function defined by A = (3/2. However, the 
correlation functions defined by the parameter A are re- 
lated one to each other by simple identities in the Fourier 
space. 

Let us consider the associated spectral functions, which 
are defined by the Fourier transforms 
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With the help of the identity 
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which is valid for all complex (3 C with Re(/3 C ) > 0, one 
computes 
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e -(0/2-\)E-(l3/2+\)E' | ( E \ \E')f 
1 



2tt 



M-u+(E-E')/h] dt 



dEdE'. 



Simple manipulations by means of the Fourier represen- 
tation of the delta function show that 

C ,aH = fie-OVa-AHi f e -^ E \(E + uh\0\E)\ 2 dE . 
Jr 

(4) 
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Therefore, the Fourier transforms of the correlation func- 
tions are non- negative distributions (in fact, they are 
thermal averages of certain power spectra). If we de- 
note the special value for A = by Go(oj), then we have 
the relationship 

Co,\(w) = e Awfi G (w), (5) 

which shows that all correlation functions defined by 
Eq. Q carry essentially the same information. It is 
not difficult to see that the spectral function Go(w) 
is symmetrical about the origin and we shall refer to 
the corresponding correlation function as the thermally- 
symmetrized correlation function. 

Berne and Harp 22 have pointed out that the computa- 
tion of thermally-symmetrized quantum correlation func- 
tions 

G (t)=tr(e-^ H 0h-^ H 0), (6) 

with C = 0/2 - it/% and % = 0/2 + it/%, might 
be an easier computational task. Certain quantities of 
physical interest, such as rates of reaction or diffusion 
constants, only depend on the value at the origin of 
the spectral functions, value that is independent of the 
particular choice of A. Miller, Schwartz, and TrompS 
have utilized this independence to point out that the 
thermally-symmetrized flux-flux correlation function has 
better mathematical properties than the Yamamoto flux- 
flux correlation function, 24 which corresponds to an av- 
erage over A on the interval (—0/2, 0/2). More recently, 
Predescu and Miller— have argued that the thermally- 
symmetrized spectral function is the one for which the 
moments (and, in general, any short-time information) 
are the most sensitive with respect to changes in the val- 
ues of the spectral function near the origin. In other 
words, for this particular choice of correlation function, 
the values of the spectral function near the origin are ex- 
pected to have the best continuity properties with respect 
to variations in the moments. 

From Eq. Q and the inverse Fourier transform for 
Eq. it follows that 

G (it) = C O -t/h(0) = / C -t/h(u)duj, 
Jr 

an equality that holds provided that t 6 (— 0H/2, 0H/2). 
On the other hand, from Eq. JSJ, we obtain the Laplace 
identity 

G (it) = f e-^G (uj)diu, (7) 

JR 

the inverse of which is the thermally-symmetrized spec- 
tral function. 

Having reached this point in our presentation, we 
pause and ask whether or not Eq. J7J uniquely deter- 
mines the spectral function. The answer is affirmative 
and follows from different arguments, all of which are 



based on the fact that the spectral function is positive. 
Thus, one could follow the path of Baym and Mermin 4 
and use positivity to argue that the integrand of Eq. Q 
is absolutely integrable. In turn, this implies that the 
correlation function Go(t) must be analytic in the com- 
plex plane on the strip defined by |Im(t)| < 0K/2. Stan- 
dard results of complex analysis then show that Go(t) 
for real t and, therefore, Go(w) are uniquely determined 
by the values of Go (it) on the interval (— 0H/2, 0H/2). 
The absolute integrability of e" flJ Go(w) plays an impor- 
tant role in the proof of uniqueness and it should not be 
taken easy. Indeed, if the integral in Eq. Q is only re- 
quired to converge in the Cauchy principal value sense, 
then there exist an infinity of solutions, of which only 
one is positive [non-positive examples are furnished by 
the Fourier transforms of Eq. The issue is rele- 

vant because both absolute integrability and analyticity 
are constraints on the set of admissible solutions that are 
very difficult to implement on a computer. By compari- 
son, enforcement of positivity is a more achievable goal. 

If we also use the a priori information that Go (it) and 
Go(to) are symmetric, the problem that must by solved 
in the context of the inverse Laplace transform method 
is: find the positive and symmetric distribution Go(u) 
that satisfies the equation 

G (it) = / cosh(ujt)G (uj)duj, (8) 

JR 

for all t G [0, 0h/2). The input data for this problem 
is usually a finite sequence of values of the imaginary- 
time correlation function on an equally-spaced grid 
{tn,j — j0h/(2n) : n > 1, < j < n}. Upon discretiza- 
tion, the functional equation exhibits multiple solutions 
and becomes determinate only upon the specification of 
an inversion algorithm. The main problem a computa- 
tional physicist has to face is that the original functional 
equation is ill-posed in the sense of Hadamard. Although 
the problem has a unique solution, this solution lacks 
continuity with the input data for virtually any compu- 
tationally reasonable topology. For example, there are se- 
quences of functions f e (t) with | f e (t) — Go (it)\/\ Go (it)\ < 
e for all t € [O,0H/2), such that the problem 

f e (t) = [ cosh(ut)G (u)dw, (9) 
Jr 

has no solutions for any e. Thus, just by mere control 
of the relative errors in the input data, we are not even 
guaranteed an inversion algorithm, much less a sequence 
of approximations to the spectral function that converges 
to Go(u) as e — > 0. 

Another approach to proving the uniqueness of the so- 
lution of the inverse Laplace transform is via moments. 17 
First, one utilizes the positivity of the spectral function 
to demonstrate that the Taylor series of the correlation 
function Go(t) about the origin has a convergence ra- 
dius equal to or larger than 0K/2. The sequence of even 
derivatives of the imaginary-time correlation function at 
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the origin reads 



D 



2 k 



d Ik G (it) 



dt 2k 



t=o 



= I G o {u)L0 2k dw 



(10) 



The odd moments are zero, by the symmetry of the spec- 
tral function. The ensuing symmetric Hamburger mo- 
ment problem is then demonstrated to be uniquely de- 
termined, thus both proving the uniqueness of the recon- 
structed spectral function and suggesting an alternative 
computational approach. Unfortunately, the inverse mo- 
ment problem also lacks continuity with the input data. 
Thus, just by controlling the relative errors for the mo- 
ments, we are not guaranteed that an inversion algorithm 
exists. Moreover, the finite moment problem may also be 
undeterminate. If determinate, the finite moment prob- 
lem is exponentially unstable. These stability issues will 
be addressed in the following section. 

We conclude this section by emphasizing again the 
requirement of positivity for the reconstructed spectral 
function. For any spectral function Go(w), there are 
modifications that satisfy both the full moment prob- 
lem given by Eq. IjlUfl and the Laplace equation given by 
Eq. (JSJ , precisely for the same input data as the physical 
spectral function. Such modifications can be obtained 
by adding some integrable and infinitely differentiable 
function that vanishes within the interval (— f3h/2, j3h/2) 
to the correlation function and then taking the Fourier 
transform. A specific example of a function that satisfies 
both the full moment problem and the Laplace equation 
is provided by the Fourier transform of 



G^\t)^Go(t) 




exp 



i-2\t\/(hp) 



if \t\ < 0h/2, 
otherwise, 



(11) 

for any arbitrary and positive parameter a. Such a 
Fourier transform always has a non-vanishing negative 
part. Of course, in agreement with Baym and Mermin's 
analyticity argument, the modification to the correlation 
function expressed by Eq. l|llf) is not analytical. Nev- 
ertheless, there is no effective procedure to enforce an- 
alyticity numerically and the positivity of the spectral 
function comes in handy. We stress that this positivity 
must be enforced to machine accuracy: for many modifi- 
cations, the negative part of the modified spectral func- 
tion appears in the high frequency region and can be a 
very small number, difficult to recognize on a plot, even if 
the correlation functions are "obviously" different. This 
is just another manifestation of the instability of the in- 
verse problems. 



III. STABILITY OF THE INVERSE FINITE 
MOMENT PROBLEM AND VERIFICATION OF 
THE MOMENT DATA 

We begin this section with a short review of the Ham- 
burger moment problem. All mathematical information 



contained in the present section can be found in stan- 
dard references on the moment problem^ A sequence of 
numbers Do, D\, . . . is called a moment sequence if there 
exists a non-negative distribution, say Go(w), such that 



D k = 



Go(uj)u k dLo. 



(12) 



Quite clearly, not all sequences of numbers are moment 
sequences. For example, any moment of even order must 
be a non-negative number. Even more, by the non- 
negativity of the distribution Go(w), we also have 

Go(oj)(^2ajLjA doj > (13) 

for all sets of number ao, a%, . . . , a n . In matrix language, 
the last inequality is equivalent to the condition that the 
Gram matrices 



E D 



A' = 



/D D 1 
D x D 2 
L>2 D« 



D 2 
D 3 



D,, 
D h 
D„ 

Do 



\ 



(14) 



are positive semi-definite (that is, their lowest eigenvalues 
must be greater or equal to zero). A standard result from 
matrix analysis says that the Hermitian matrix A' n is 
positive semi-definite if and only if all determinants |AJJ 
with < k < n are non-negative. Therefore, a necessary 
condition for a sequence of numbers to be a moment se- 
quence is that the matrices A' n are positive semi-definite 
for all n > or that the determinants \A' n \ are non- 
negative for all n > 0. Hamburger has demonstrated 
that these conditions are also sufficient for a sequence 
of numbers to be a moment sequence. In addition, he 
has shown that, given a finite sequence Do, D\, . . . , Z?2 n , 
the positive semi-definiteness of is sufficient for the 
finite moment problem to have at least a solution (obvi- 
ously, such a solution is rarely unique). If the quantities 
Di, Z?3, . . . , Z?2n-i are zero, then there exists at least one 
symmetric solution. 

As shown by Eq. flU , the Gram matrices A' n have a 
very special structure: the skew-diagonals are made up 
from identical elements. Such matrices, whether Gram or 
not, are called Hankel matrices. For the symmetric Ham- 
burger moment problem, the skew-diagonals correspond- 
ing to moments of odd order are zero. The importance 
of the positive semi-definiteness of the matrices can 
be understood in the context of Theorem 3 of Ref. ITtI 
which is a standard convergence theorem in probability 
theory. Namely, assuming that we are given a collection 
of moment sequences D Q n \ D^\ . . . , D g„ (the low-rank 
terms of which are allowed to change with n for gener- 
ality) and letting Gg'(u) and G^^t) denote the associ- 
ated spectral and correlation functions, respectively, the 
convergence 



lim D 



(n) _ 



D k , V k > 



5 



implies 



lim G { o\t) =G (t), VtE 



This result appears to contradict our previous assertion 
that the moment problem lacks continuity with the in- 
put data. To the contrary, the theorem provides a means 
of approximating the exact correlation function. The 
explanation is that Theorem 3 requires the input data 
Dq^ , D± , . . . , £^2„ to be finite moment sequences and 
it is this requirement that lacks continuity with the in- 
put data. More exactly, for any e > 0, there exist 



a rank n and numbers , D 



(n) n (n) 



l-Dfe" 5 - Ac I < V < k < 2n, yet the new data 

are not a moment sequence (their Gram matrix is not 
positive semi-definite). Let us consider a particular case 
where e = 0.01. Thus, we know all the moments with 1% 
relative accuracy. Is this enough to be able to generate a 
good approximation to the correlation function? The an- 
swer is no. As the results in the remainder of the present 
section show, it is very likely that the data we posses do 
not form a moment sequence, even for moderately large 
n. 

For our symmetric problem, the moments of odd order 
are zero and, therefore, their value is exactly known. In 
agreement with the hypothesis of Theorem 3, we require 
of any computational procedure to be able to provide 
the even-order moments with controlled relative error. 
By making this relative error small, we may assume that 
the even-order moments are positive. According to the 
discussion in the preceding paragraph, for the reconstruc- 
tion algorithm to converge to the exact result in the limit 
that the relative error for the even-order moments con- 
verges to zero, it is sufficient that the inequalities 



,D%£ such that 



ajD j+k a k 




(15) 



are satisfied for all numbers do, ffli, ■ ■ ■ , a n . By the posi- 
tivity of the quantities D 2 j, this condition is, of course, 
equivalent to the one provided by Eq. I|13|) . However, 
it also takes into account the fact that the relative er- 
rors of the moments D 2 i are controlled. By making the 
substitution aj = a'j/ \J~th~j in Eq. Ijl5|l . we see that the 
above inequality is equivalent to the condition that the 
Hermitian matrices A n of entries 

(A„)j,fe = Dj +k J \jD 2j D 2k , < j, k < n 

are positive semi-definite. We summarize the findings 
obtained so far in the present section in the following 
theorem, which gives sufficient criteria for the existence 
of well-posed inversion algorithms. 



(«) n(") 



Theorem 1 For each n> 1, let Dq ,D 

a finite sequence of positive moment data. Let 



be 



(n) 



0. 



for k = 1, 3, . . . , 2n — 1, and assume that the Hermitian 
matrices A„ of entries 



(A n ) jtk = D 



sa/( 



D v ' D 



2k 



1/2 



0<j,k<n (16) 



are positive semi- definite. Then there exists at least one 
symmetric trial spectral function Gq^ (lu) of even mo- 
ments D { n) , D ( 2 n) Ain • Moreover, with G^ (t) de- 
noting the associated trial autocorrelation function, the 



convergence 



implies 



lim D 



lim G { o\t) = G (t), Vie 



(«) _ 



D k , V k > 



(17) 



(18) 



The upper index (n), which was needed in the formula- 
tion of the theorem, will be dropped from now on. We 
shall use the notation D 2 k for the moment data and un- 
derstand that they are subject to both systematic and 
statistical errors. 

In view of the above theorem, it is quite unfortunate 
that the matrices A„ are ill-conditioned (although they 
are better behaved than the matrices A^). In a research 
born out of frustration with the numerical instabilities of 
an otherwise reasonable algorithm, Tyrtyshnikov 27 has 
demonstrated that the condition number (the ratio be- 
tween the largest and the smallest eigenvalues) of any 
positive semi-definite Hankel matrix grows at least ex- 
ponentially. As adapted to our problem, Tyrtyshnikov's 
result states that 



k(a;)>3-2»- 



(19) 



A tighter bound has been given more recently by 
BeckermannpS who has demonstrated that 



k (a;j > 7 "~7(i6n), 



(20) 



for n > 3. The quantity 70 « 3.210 is related to the 
so-called Catalan series, but the exact value is not im- 
portant for our purposes. Nevertheless, Beckermann has 
demonstrated that this is the best estimate for the min- 
imal value of the exponential factor 70. Thus, there are 
positive semi-definite Hankel matrices for which the ex- 
ponential growth is exactly Jq . For most other applica- 
tions, the exponential factor 



7 = lim [ K (A' n )j 



l/r, 



(21) 



is larger than 70 ~ 3.210 and may itself increase to infin- 
ity. 

Although the matrices defined by Eq. (|l(j[l are not Han- 
kel, the extreme ill-conditioning of the matrices A^ car- 
ries over to the matrices A„. A simple example will con- 
vince the reader of this. Consider the following spectral 
functionM 



G F (u>) 



1 Mft5 



(3h 2tt 



-K 



(22) 



6 




FIG. 1: Asymptotic behavior of the quantities k(A„) 1 ^" for 
the flux-flux spectral function of a free particle. The asymp- 
totic behavior demonstrates the exponential instability of the 
matrices A„. The condition number of the stability matrices 
increases roughly as 2.3™ for large n. 



where K\ (x) denotes the respective modified Bessel func- 
tion of the second kind. This is the spectral function asso- 
ciated with the thermally-symmetrized flux-flux correla- 
tion function (the so-called Miller, Schwartz, and Tromp 
correlation function^ or MST for short) for a free par- 
ticle 



G F (t) 



1 



Ph [ t 2 + (/3ft/2)2] 3 / 2 ' 



(23) 



Formal differentiation of the correlation function at the 
origin shows that the even-order moments of the MST 
spectral function for a free particle are given by the equa- 
tion 



1 2k + 1 (2fc)P 
D2k = — - 



(24) 



7r(/3ft) 2fe + 2 kl 2 ■ 
Using the above formula, one can set up the Hcrmitian 



matrices A„, diagonalize them, and compute their con- 
dition number n(A n ). As apparent from Fig. ^ the 
quantity ^(An) 1 /" converges to a constant, the estimated 
value of which is 7 ~ 2.3. This implies that the sequence 
of matrices A„ is exponentially unstable. As Eq. (0J 
with A = suggests, the tail of the spectral function de- 
cays like e ~0 h ^/ 2 (as exponential order) for all thermally- 
symmetrized spectral functions. Because the relative val- 
ues of the high-order moments depend only on the prop- 
erties of the tail of the spectral function, we suggest that 
the asymptotic value of 7 « 2.3 may be characteristic 
of all spectral functions. The results computed for the 
symmetric Eckart barrier (they are presented in Fig. 0) 
seem to support the suggestion. 

Why is the exponential instability of the matrices A„ 
so important? It shows that, just by controlling the rel- 
ative errors of the even-order moments (that is, just by 
controlling the relative errors of the entries of the ma- 
trices A n ), it is very likely that we cannot ensure the 
positive semi-definiteness of the matrices A„, even for 
moderately large n. As it is well known, the condition 
number of a matrix controls the relative errors in the 
values of the determinants |A n | with respect to the rela- 
tive errors in the entries of the matrix. Roughly speak- 
ing, log 10 (K(A„)) represents the number of exact digits 
with which the entries of the matrix A„ must be known 
in order to guarantee that the determinant |A n | is still 
positive. Therefore, the positive semi-definiteness of the 
matrices A n must be a "built in" feature of the compu- 
tational procedure. 



IV. MOMENT ESTIMATING FUNCTIONALS 

As discussed in Ref . Il7l the moments can be computed 
by Monte Carlo integration as averages of some estimat- 
ing functionals. A typical average is expressed by the 
equation 



1 d k 



N F Nf dt k 



G F (it) 



f s dxdzEE'e- ||z||2 e~ (/3/2) [-'« 1 V{*+°a™+a a Bl)du+Q v(x+a a zu+* B°')du] d^^ ^ Z) B°,B°' 
f dx<feEE'e-|HI 2 e~ (/3 < /2) [-'° 1 V (^+^u+a B°)du+J^ V(x+a zu+<T B°')du] 



(25) 



The significance of the various terms can be found in 
the cited reference and will be partially explained below. 
Nonetheless, Eq. I|25|l is sufficient to point out one of the 
main numerical difficulties of the present approach: due 
to the extreme complexity of the estimating functional 



Tt (x, z, B®, B®'j , the differentiation against the param- 
eter t cannot be done analytically. As already mentioned 
in the introduction, the differentiation can be performed 
by finite difference. Such an approach, however, requires 
high precision evaluation of the potential function enter- 
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ing the expression of the estimating functional. 

In the mathematical literature, it is well known that 
a superior technique for performing numerical differenti- 
ation is provided by the so-called pseudo-spectral meth- 
ods. The functional Tt (x, z, B®, B°'J is infinitely differ- 
entiate on the interval t G (— (3H/2, /3H/2) provided that 
the potential function is also infinitely diffcrcntiable. The 
functional can be approximated on some compact inter- 
val [-9/3h/2,9/3h/2}, < 9 < 1, by Chebyshev polyno- 
mial interpolation (or, following certain transformations, 
one can also perform a trigonometric interpolation) . The 
resulting Chebyshev polynomial can then be differenti- 
ated analytically. It can be demonstrated that, when 
differentiating such an interpolating polynomial, the er- 
ror committed decays to zero exponentially fast for func- 
tional that are analytical in t (that is, for analytical 
potentials) and faster than any polynomial for infinitely 
differentiable functionals (potentials). This should be 
compared with the polynomial decay obtained by finite 
difference. A short exposition of these results and some 
elegant proofs can be found in Ref. |3(1 

The choice of the parameter 9 is crucial in obtaining ac- 
curate estimates of the derivatives. If 9 is too small, then 
the input information represented by the values of the es- 
timating functional Tt fx, z, £?°, B+'^j for the Chebyshev 

knots in the interval [—9/3H/2, 9(3h/2] are highly redun- 
dant and the necessary precision is very high. Therefore, 
the parameter 6 should be chosen large, as close as pos- 
sible to the value 1 (say 9 = 3/4). Even so, the number 
of Chebyshev coefficients that are accurately computed 
is not very great and depends on the precision of the ma- 
chine. As a rule of thumb, one can rely upon 8 to 16 coef- 
ficients in single precision and 16 to 32 coefficients in dou- 
ble precision. As recommended in Numerical Recipes,^- 
a technique that improves the quality of the polynomial 
interpolation is to truncate a higher-order interpolating 
polynomial to a lesser degree. We shall refer to such an 
interpolation polynomial as a regressed polynomial. 

Given the limitation due to the finite precision with 
which the potential can be evaluated, the pseudo-spectral 
technique may fail to provide adequate estimates for the 
derivatives if the estimating functional and its derivatives 
are not easily approximated by a low-degree interpolant. 
Unfortunately, this is the case for the present computa- 
tional task. The culprit is the Boltzmann factor that en- 
ters the definition of the estimating functional. To under- 
stand how this factor enters our equations, we review the 
definition of the flux-flux estimating functional^ The 
surface through which the flux is computed is assumed 
to be a plane of equation xi = 0. Thus, the space iS ap- 
pearing in Eq. I|25ll is the subspace of K d x M. d of equation 
xi = and %\ = 0. The quantities B® and B® are in- 
dependent d-dimensional standard Brownian bridges (d- 
dimensional vector valued stochastic processes, the com- 
ponents of which are independent one-dimensional stan- 
dard Brownian bridges). We also define the entities f} t , 

1 /2 

a t , and a± t as /3 t = j3/2 + t/H, a t = (h 2 p t /m ) 



and a±t — <Jt<T-t/&o, respectively. Finally, we let V(x), 
V'(x), and V"(x) represent the potential and its first 
and second order partial derivatives against the reaction 
coordinate Xi. 

The following notation, additional to what has been in- 
troduced in Ref. 0, defines several action-like variables 
which will constitute the basic entities that we interpo- 
late: 

S t (x,z,B°) =fa [ V{x + a ±t zu + a t B Q u )du, (26) 
Jo 



S't a {x, z, B°) = fit / V (x + a± t zu + o t Bl) udu, 



(27) 



S' t ' b (x,z,.B°) =(3 t f V (x + a± t zu + a t B°) (l-u)du, 
Jo 

(28) 



and 



S t "(x, z , B°) = ft [ V" (x + a± t zu + a t B°) u(l - u)du. 
Jo 

In terms of the action-like variables, we define 



V / a_t at 



S' t > a (x,z,B°:)-S'l a t (x,z,B°) 
^(x,z,Bf)-^(x,z,J]!)] (30) 
S' t '( X ,z,B°J)-S': t (x,z,B°) 



and 



A t (x,z,5°) = - [S (x,z,B°)-S t (x,z,B°)] . (31) 
Then, as shown in Ref. 0, 

^(x, Z ,^^)^^{^ t (x,z,^,^') 
xeW2 )[A t (x, z X)+A- t (x,,, B r)] + jrO (x !Z)jB )jB 0^ (32) 
xe (/3/2)[A_ t (x,z,B°)+A t (x,z,B°')]| _ 

The weighting factors of the type 

e (0/2)[A t (x,z,B°)+A- t (-K,z,B°')] 

may vary quite violently for low enough temperatures 
and cannot be readily approximated with a low-degree in- 
terpolant even for smooth potentials. On the other hand, 
the action-like quantities do not vary violently with t even 
for low temperatures. The extent to which they vary is 
controlled by the values of the potential or its derivatives 
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(the smoothness of these functions is assumed). It is, of 
course, these quantities that we intend to approximate 
by pseudo-spectral methods. 

For this purpose and in order to take advantage of the 
whole array of values t £ [— /3H/2, /3K/2], we make the 
substitution 



t 



= sin 



-(1 



(33) 



and regard the action-like entities as functions of <f> on 
the interval [—1, 1]. With the notation a — {h 2 P/m®) 1 / 2 , 
the equalities 



0t 


= /3sin 




= a sin 




1 . 




= —asm 
2 



4^ 
t(1 



(34) 



show, for example, that the action St can be regarded as 
the function of <j> given by the formula 



^(x )Z ,S°)=/3sm [-(1 



r ( x + 2 azu 



x sm 



-(1 
2 V 



cjB° sin 



-(1 

4 V 



du. 



This function is infinitely differentiable on the interval 
[—1,1] whenever the potential is so. Therefore, the 
action-like functionals can be Chebyshev approximated 
on this interval faster than any polynomial. An alterna- 
tive interpolation procedure utilizes trigonometric poly- 
nomials and is based on the observation that the action- 
like variables are periodical in (f>. Unfortunately, the pe- 
riod is 8 and the functionals must be sampled on the 
larger interval <fi G [—4,4], which is more costly, because 
one needs roughly 4 times more points. Although the 
trigonometric (Fourier) interpolation has the nice prop- 
erty that it becomes exact for potentials that are poly- 
nomials, we did not notice any significant advantage over 
the Chebyshev interpolation in more realistic numerical 
tests (the two techniques behave in a similar fashion for 
comparable meshes of the interpolatory knots). 

For some order of approximation n, let <fik — cos(7r(fc — 
0.5)/n), k — 1,2, ... ,n, be the nodes of the Chebyshev 
polynomial T n {x). For j = 0, 1, . . . , n— 1, the Chebyshev 
coefficients are given by 



c,(x, z, B°) ^-^S^ (x, z, B°) cos 



k=l 



nj(k - 0.5) 



and we have the following approximation formula 



(35) 



. n— 1 

S (x,z,£?°) « - C o(x,z,B°) + ^ Cfc (x,z,i?°)T fc (0). 



k=l 



(36) 

The right-hand side expression in Eq. (|36|l is called the 
Chebyshev interpolation polynomial of rank n — 1. This 



polynomial is the unique polynomial of rank n — 1 that 
coincides with ^(x, z, B®) for all n interpolatory knots 

Of course, interpolation polynomials must be com- 
puted for all action- like functionals described by Eqs. I|26|) 
to l|29|) . The computation of the Chebyshev coefficients 
can be performed in a fast and stable way by cosine 
fast Fourier transform (FFT)rfJi We recommend the use 
of such a transform not so much for reasons of com- 
putational efficiency as for reasons of accuracy. If e is 
the floating-point relative precision of the machine, the 
relative error for the Cooley-Tukey FFT algorithm is 
0(elog(n)), compared to 0(en 3 ' 2 ) for the direct matrix 
multiplication technique^ The numerical evaluation of 
the action-like functionals from their Chebyshev coeffi- 
cients is to be performed by the Clenshaw recurrence 
formula, which is documented in Numerical Recipes^i 

Replacing the action-like functionals with their Cheby- 
shev approximation in Eqs. 13U|) to (|32|) . we obtain a 

non-linear approximation T^, ^x, z, B®, B®'^j of the esti- 
mating functional in terms of the variable (f>. In terms of 
the variable t, one has 



T t (x, z, Bl B°J) w P m (x, z, Bl B°J 



(37) 



where 



= -i 



7T 

Ai 



In 



t 

hp 



t 

hp 

+ i 



t 

hp 



(38) 



is, of course, the appropriate solution of Eq. H33JI . 
Eq. Q38JI. which involves the use of complex numbers, 
already betrays the approach we shall utilize to compute 
the derivatives of the estimating function at the origin: 
contour integration of the complex extension of the non- 
linear Chebyshev approximation in the complex plane. 
The numerical algorithm utilized is due to Lyness^ and 
is summarized in the following paragraph. 
The Cauchy integral formula 



k\ 

27ri 



B°J 



dt 



(39) 



provides a way to compute the derivatives of an analytical 
function by computing integrals. The contour C must be 
a closed curve that contains the origin in its interior. It 
will be taken to be a circle of radius r E (0, Ph/2) centered 
about the origin. Eq. H39|) becomes 



-2-Kik8 



T 



U^Bl,Bl' 



dO 



(40) 



and shows that the problem of computing derivatives is 
equivalent with that of evaluating the Fourier coefficients 
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of a periodic function. The one-dimensional integral in 
Eq. 1)41) [) is to be computed by trapezoidal quadrature. 
We have the approximation 

" T^E^ tti/m W (x,.,a°,fl2') ■ (41) 
r K m * — ' V / 

where the summation is best executed by FFT. Lyness 33 
has demonstrated that this simple-looking algorithm is 
numerically stable and converges exponentially fast with 
respect to the number of quadrature knots m. For this 
reason, the algorithm is known as Lyness' method. Al- 
though criteria for choosing optimal values for the radius 
r are known, numerical tests show that the ad-hoc value 
of r = /3H/4 produces excellent results. 

We conclude the present section by pointing out an 
issue of convergence and precision. Unless the potential 
function is analytical, one cannot extend the estimating 

functional Tt fx, z, B®, B®'^j to the entire complex plane. 
Thus, the non-linear Chebyshev approximation, although 
convergent on the real interval (— /3H/2, (3H/2), diverges 
on the complex disk of radius 0H/2, as the rank n of the 
interpolating polynomial goes to infinity. We take advan- 
tage of the rapid convergence and the stability of Lyness' 
method to compensate for this divergence. However, to 
also compensate for the loss in precision in the evaluation 
of the integrand, one may need to evaluate the Cheby- 
shev polynomials as well as the final non-linear approx- 
imation in higher precision. Numerical tests show that 
this divergence is very weak and that the need for higher 
precision seems only theoretical. With the increase in the 
rank of the Chebyshev approximation beyond a certain 
level, we do notice a sudden divergence, but this is caused 
by the limited precision with which the high-order coef- 
ficients are determined and appears even for analytical 
potentials. The convergence of the Chebyshev approxi- 
mation is normally so fast that it rapidly feels the lack of 
smoothness of the action due to the finite precision of the 
machine. Thus, the high-order coefficients, instead of de- 
caying steadily to zero, remain roughly constant as mag- 
nitude. One counteracts this artifact by truncating the 



Chebyshev series to a safe number of coefficients (usually 
less than 16 in single precision and less than 32 in double 
precision). The sudden divergence can be easily spot- 
ted by comparing the values of Do computed with two 
slightly different estimators: the first one utilizes the es- 
timating functional Tq (x, z, B®, B®'^j (more precisely, to 

account for the case when the polynomials are regressed, 
the value at the origin of its non-linear approximation 
in terms of Chebyshev polynomials), whereas the second 
one utilizes the Cauchy contour integral 

^ m 

-YsFre^l™ (x,Z,D?,D?') . 

i=i 

The agreement for the computed values of Do must be 
better than the minimal accuracy that would guarantee 
the positive semi-definiteness of the stability matrices. 



V. A NUMERICAL EXAMPLE: THE 
SYMMETRIC ECKART BARRIER 

In order to demonstrate the capabilities of the present 
technique, we compute the first 20 even-order moments 
of the flux-flux spectral function for a symmetric Eckart 
barrier at the low temperature of T = 100 K. The pa- 
rameters for the Eckart barrier are those also utilized in 
Ref. ^3- The potential is 

V(x) = V sech(ax) 2 , (42) 

with Vq = 0.425 eV, a = 1.36 a.u., and tuq = 1060 a.u. 
The parameters for the barrier are chosen to correspond 
approximately to the H + H2 reaction. As discussed in 
Section III, the necessary and sufficient criteria that guar- 
antee the existence of a convergent reconstruction algo- 
rithm are a) positivity and control of the relative errors 
for the moments and b) positive semi-definiteness of the 
matrices A„. 

Let us summarily describe the main features of the 
Monte Carlo path integral technique utilized. According 
to Eq. H25|) . what we actually compute are the ratios 
D 2 fe/A/>, where 



Af F = 



1 



1 



d-i 



Sttjtiq \2tT(Tq / 



(43) 



Although the normalization coefficient J\fp can also be 
determined by Monte Carlo integration, its value is irrel- 
evant for the present study. The path integral technique 
utilized is based on the fourth-order short-time approxi- 
mation introduced in Ref. |3 A Trotter index of 16 has 



been employed, for a total of 64 path variables. In or- 
der to diminish the amount of correlation in the Monte 
Carlo chain, we have utilized the so-called fast sampling 
algorithm considered in Ref. [3^ as well as the parallel 
tempering technique. The parallel tempering exchanges 
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have been performed with replicas of temperatures ar- 
ranged in geometric progression up to 5000 K. A total 
of 10 million complete sweeps through the space of path 
variables have been made. These sweeps have been di- 
vided in 100 blocks. 

Naturally, the Monte Carlo simulation consists of two 
parts: sampling and accumulation of averages for the 
estimating functionals. Due to the nature of the fast 
sampling algorithm, which organizes the path variables 
in 2 + log 2 (64) = 8 layers with the variables from the 
same layer sampled separately (and independently), the 
computational cost for a sweep is 64 • 8 calls to the poten- 
tial function. The computational cost for the estimating 
functionals is 64-3-32 calls. The factor 3 comes from the 
three different types of functions that are called [V(x), 
V'(x), and V^"(x)] whereas the factor 32 represents the 
number of Chebyshev knots. We mention that even for 
such a simple analytical potential, the largest part of the 
computation is spent on the evaluation of the potential 
and its derivatives and not on performing the numer- 
ical manipulations considered in the preceding section. 
In order to dedicate comparable amounts of time to the 
sampling and estimator evaluation parts, we have eval- 
uated the estimators every 3-32/8 = 12 sweeps. We 
made this analysis just in order to warn the reader about 
the magnitude of the computational cost incurred by the 
evaluation of the estimating functionals for derivatives. 
Thus, it is not worth to accumulate averages after each 
sweep, especially given the usually large correlation of the 
Monte Carlo sampler (it is very rare that the correlation 
times are smaller than 12 in realistic simulations). 

As mentioned before, the main purpose of the present 
development is to provide a technique that is capable of 
producing accurate estimates for the high order deriva- 
tives without the need to evaluate the potential function 
in arbitrary precision. Therefore, the potential function 
and its derivatives utilized in the construction of the es- 
timating functionals have been evaluated in single pre- 
cision (about 7 significant decimal digits) only. This is 
also the precision with which the Chebyshev coefficients 
are evaluated. On the other hand, the evaluation of the 
Chebyshev polynomials, the contour integration in the 
complex plane, and the accumulation of averages have 
been performed in quadruple precision (about 31 signifi- 
cant decimal digits). The sampling part of the simulation 
has been conducted in double precision (about 15 signif- 
icant decimal digits). To complete the description of the 
algorithm, we mention that the Chebyshev polynomials 
have been regressed to the first 16 coefficients and the 
contour integration has been performed in 64 quadrature 
knots. 

We accumulate the averages in quadruple precision in 
anticipation of the loss of precision due to the exponential 
instability of the matrices A„. We stress that, although 
the moments Z?2fc are determined with a precision of a few 
digits only, an important amount of information resides 
in the remaining imprecise digits. To understand how 
this may happen, let us analyze the problem of comput- 



ing the moments if the spectral function is analytically 
known and can itself be sampled. Thus, we compute av- 
erages of the type 

l x k = D u /D = J G F {Lj)uj k duj I J G F (u)duj (44) 

for k = 0, 1, ... by Monte Carlo integration. The actual 
ratios in N sample points are 

1 N r 
= iE^= / PN(^ k dcj, (45) 

where 

1 N 

Pn(uj) = —^2S(uj -ujj) (46) 
i=i 

is some discrete measure that depends upon the history 
of the Monte Carlo integration. We notice that irrespec- 
tive of what this discrete measure is, the sequence of 
numbers {fj, k ;k — 0, 1, . . .} is a moment sequence be- 
cause it comes from some positive measure. It therefore 
satisfies the requirements of positive semi-definiteness of 
the stability matrices A„ regardless of the actual num- 
ber of samples N. It also exhibits the same instability 
issues as the original problem. Theoretically, if we deter- 
mine the whole sequence of moments {(-i k ; k = 0, 1, . . .} 
for some fixed N with arbitrary precision, we can recon- 
struct the measure that has generated the sequence, that 
is, the distribution pj^{ui). But if these statistically in- 
exact moments are not known with sufficient precision, 
the exponential instability will cause them to loose the 
crucial property that they form a moment sequence. 

The lesson we learn from the preceding exposition is 
that we cannot compute moments of different orders in 
independent Monte Carlo runs, because it is difficult to 
attain the precision necessary to ensure the positive semi- 
definitcness of the stability matrices A„. In this paper, 
we declare ourselves content with the simultaneous com- 
putation of the moments in the same Monte Carlo simu- 
lation. We mention, however, that this is not necessarily 
without penalty. If, for example, one is interested in 
evaluating the tail of the spectral function [which decays 
like e - ' 3 '""'/ 2 , as Eq. J3J with A = suggests], then there 
are going to be problems related to the fact that the tail 
of the distribution is infrequently sampled by the Monte 
Carlo walker because of its exponentially vanishing sta- 
tistical weight. Accordingly, the information contained 
in the moments determined by Monte Carlo integration 
poorly reproduces the properties of the tail. Thus, our 
assumption is that we arc interested in the properties 
of the spectral function near the origin or in regions of 
important statistical weight. 

Under this assumption, the poor statistics for the com- 
putation of high order moments mentioned in the pre- 
ceding paragraph is harmless. It also serves to show that 
the relationship between the quality of the reconstructed 
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TABLE I: Moments of even order and their percent relative statistical errors determined by Monte Carlo integration. A number 
of log 10 [fi(A2o)] ~ 8 significant figures are necessary in order to prevent the loss of the positive semi-definiteness property of 
the stability matrices. We give the results with one significant digit more than the minimal requirement so that the reader may 
verify that the matrix A20 is positive definite. 



k 





1 


2 


3 


4 


5 


6 


D 2k /J^F 


5.43598845 • 10 1 


2.25127499 • 10~ 4 


3.79457212 • 10~° 


1.31959326 ■ 10~ 13 


7.49921061 ■ 10~ 18 


6.09205035 • 10~ 22 


6.48529540 • 10~ 26 


Error 


0.5% 


0.7% 


1.2% 


2.4% 


4.3% 


6.7% 


9.5% 


k 


7 


8 


9 


10 


11 


12 


13 




8.53623493 ■ 10~ 30 


1.33546663 • 10~ 33 


2.41660791 • 10~ 37 


4.96421869 ■ 10~ 41 


1.14347963 ■ 10~ 44 


2.93267228 • 10 -48 


8.35072926 • 10 -52 


Error 


12.6% 


15.7% 


18.8% 


21.5% 


23.7% 


25.2% 


26.0% 


k 


14 


15 


16 


17 


18 


19 


20 


D 2k /^F 


2.64105405 ■ 10~ 55 


9.29940331 • 10~ 59 


3.65639217 • 10~ 62 


1.60895318 ■ 10~ 65 


7.92788909 ■ 10~ 69 


4.36750842 • 10~ 72 


2.68123444 • 10~ 75 


Error 


26.2% 


25.9% 


25.4% 


24.1% 


24.3% 


23.8% 


23.5% 



spectral functions and the relative errors of the moments 
is far from being a linear one. However, we do need 
to worry about the fact that the estimating function- 
als for the high order derivatives appearing in Eq. I|25() 
are generally not positive quantities. As such, there is 
the real possibility that, due to large cancellations be- 
tween the positive and negative parts, the computed se- 
quence Do, D\, Dii ■ ■ ■ may not be a moment sequence. 
In fact, taking into account the exponential instability 
of the matrices A„, one needs to worry about the lim- 
itations caused by the utilization of a finite Trotter in- 
dex as well as about the inherent numerical limitations 
of the pseudo-spectral methods utilized. All these sys- 
tematic errors are additional to the statistical errors. It 
appears then quite surprising that the present approach 
is extremely capable in this respect. The reader may use 
the results in Table [I] to verify that the stability matrix 
A20 is, indeed, positive definite (due to their structure, 
the stability matrices A„ for n = 1, 2, ... 19 are also pos- 
itive definite). We mention that we have verified the 
positive defimteness of the matrices A20 (by performing 
a Cholesky decomposition^!) not only for the final data, 
but also for the individual averages collected for each of 
the 100 blocks in which the Monte Carlo simulation has 
been divided! Why the algorithm is so capable in deal- 
ing with the sign problem is something that the author 
cannot explain at the present time. 

As demonstrated by the results in Table U the relative 
statistical errors increase quickly with the order of the 
moments and reach a plateau at about 0.25 — 0.26. This 
behavior seems to be caused by poor statistics, in a way 
that is perhaps similar with the previously discussed case, 
where the spectral function is sampled directly. We have 
noticed that, for k > 6, the block averages fail to become 
independent and the Monte Carlo correlation times are 
comparable to the length of the simulation. It is there- 
fore of certain interest to design alternative sampling 
techniques that would improve the statistics by use of 
suitable importance functions. However, the techniques 
should preserve the property of positive semi-dcfinitencss 
of the stability matrices. Such a task appears formidable 
because the condition numbers of the stability matrices 
increase exponentially, according to the law 2.3™. This 




FIG. 2: Asymptotic behavior of the quantities K(A n ) 1/,n for 
the flux-flux spectral function of the Eckart barrier. As for 
the free particle case, we notice that the matrices A n become 
exponentially unstable. The condition number of the stability 
matrices increases roughly as 2.3" for large n. 



exponential instability is apparent from Fig. [5] 



VI. SUMMARY AND CONCLUSIONS 

We have provided an in-depth analysis of the prob- 
lem of constructing estimators for the purpose of com- 
puting moments of spectral functions by path integral 
simulations. The estimators are constructed by formal 
differentiation of a certain estimating functional against 
the imaginary time. We have argued that the numeri- 
cal differentiation can be more successfully implemented 
by means of pseudospectral methods, in a way that uti- 
lizes information from the entire interval (— (3H/2, f3h/2). 
The algorithmic detail that leads to robust numerical ap- 
proximations is the fact that the action-like functionals 
and not the actual estimating functional are interpolated. 
The derivatives at the origin can be computed from the 
ensuing non-linear approximation in a fast and stable 
way by contour integration in the complex plane, with 
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the help of Lyness' method. 

We have improved upon the convergence results of 
Ref. ^3 by stating Theorem 1, which provides the nec- 
essary and sufficient conditions for the existence of con- 
vergent reconstruction algorithms. The hypothesis of the 
new theorem is now based on assertions that are verifi- 
able solely from the computed data. In particular, we 
have pointed out that the existence of inverse algorithms 
is inherently related to the positive semi-definiteness of 
certain matrices, which, however, prove to be exponen- 
tially unstable. Although the main statements that lead 
to Theorem 1 are well-known results from the mathe- 
matical literature, we believe it is worth having them 
systematized in a single statement, for the benefit of the 
readers. 

One of the main assumptions made throughout the 
present work is that the potential function is infinitely 
differentiable. This assumption is not necessarily a lim- 
itation if one takes into account the existence of the 
partial averaging technique^ which replaces the action- 
like variables by smooth versions obtained by convolution 
with Gaussians of certain widths. Such convolutions are, 
of course, differentiable infinitely many times. 37 For most 
practical applications, the value of the Gaussian width 
remains orders of magnitude larger than the resolution 
capabilities of the working precision. In other words, 
the interpolation technique utilized still "sees" a smooth 
functional even for the largest numbers of path variables 
that make the approximation convergent for all practical 
purposes. For these reasons, the partial averaging tech- 
nique can be thought of as the natural setting for the 
implementation of the present approach. 



Several questions related to the moment approach re- 
main to be answered. The first one asks for an explana- 
tion of why the pseudo-spectral technique utilized leads 
to estimators that, upon largely inaccurate Monte Carlo 
integration, still produce a sequence of moments. By 
"largely inaccurate," we mean that the statistical errors 
are orders of magnitude larger than the working pre- 
cision required by the observed exponential instability 
of the matrices A n . Even the systematic errors intro- 
duced by the Chcbyshcv approximation are significantly 
larger than the required precision. A second question 
asks why the sign problem that should have ruined the 
Monte Carlo simulation is actually very mild. A third 
question is related to the development of sampling tech- 
niques, perhaps by means of suitable importance func- 
tions, that would alleviate the poor statistics associated 
with the computation of high-order moments, yet would 
preserve the positive semi-definiteness of the stability ma- 
trices. A final task calls for the development of actual 
reconstruction techniques of the spectral functions from 
their moments and for studies of the suitability of the 
various techniques for specific problems. 
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